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I.  INTRODUCTION 


A.  A  GENERAL-PURPOSE  FINITE  ELEMENT  PROGRAM 

A  general-purpose  finite  element  program  should  be  able  to  solve  a  variety  of 
problems  from  the  number  of  disciplines:  linear  and  non-linear,  static  and  dynamic 
problem  of  elasticity,  fluid  mechanics,  heat  transfer,  etc.  and  can  solve  problems  of 
large  size  involving  a  variety  of  elements. 

A  general  program  is  going  to  be  voluminous  and  complex.  It  is,  however, 
desirable  that: 

•  its  logic  be  easily  understood; 

•  one  or  many  of  its  parts  be  easily  modifiable; 

•  it  offers  possibilities  to  tailor  its  facilities  for  the  solution  particular  classes  of 

problems. 

The  program  should  have  a  modular  structure,  with  the  modules  made  as 
independent  from  one  another  as  is  practicable.  The  following  modular  operations  are 
recognized: 

1.  Problem  Definition  (data  base): 

-  node  coordinates  and  element  connectivities; 

-  nodal  element  properties; 

-  boundary  conditions. 

2.  Element  Computations: 

-  integration  points  and  associated  weights; 

-  interpolation  functions  and  their  derivatives; 

-  Jacobian  matrix,  inverse,  and  determinant; 

-  element  matrices  and  vectors:  [k|,  [m],  { f} ,  etc. 

3.  Assembly  Operations: 

-  assemblage  of  master  matrices  and  vectors,  |K],  |M|,  etc. 

4.  Solution: 

-  factorization  of  master  matrices  and  solution  of  equations. 

5.  Result: 

-  output  of  nodal  variables  and  other  calculated  quantities:  gradients, 
reactions,  etc. 


Subroutines  implementing  the  various  operations  described  above  are 
contained  in  all  finite  element  codes.  The  flow  of  information  between  these 
operations  is  problem  dependent;  linear,  non-linear,  static,  and  dynamic  problems  all 
require  logic  of  their  own. 

B.  MEF  PROGRAM 

A  program  of  medium  complexity,  called  MEF,  implementing  the  techniques  of 
general-purpose  program  that  can  solve  a  large  variety  of  boundary  value  problems  of 
mathematical  physics.  It  is  written  in  FORTRAN  IV  and  can  be  easily  adaptable  to 
various  computers. 

The  main  program  controls  the  flow  of  all  information  through  the  functional 
blocks  by  transferring  control  to  a  subroutine  called  BLNNNN  when  the  block  calling 
card  NNNN  is  encountered  in  the  input  file.  The  subroutine  BLNNNN  then  performs 
preliminary  functions  such  as  logical  unit  identification,  and  reading  of  control 

parameters  for  the  creation  of  various  files  and  tables.  The  subroutine  then  calls 

subroutine  EXNNNN.  In  all  cases,  subroutine  BLNNNN  provides  appropriate  default 
parameters  which  will  be  overridden  by  user  values  if  specified.  Subroutine  EXNNNN 
then  performs  the  major  operations  of  the  block  by  calling  on  the  needed  subroutines 
in  the  MEF  library.  The  above  protocol  holds  for  all  blocks  except  STOP,  COMT, 
and  1MAG.  All  the  functions  of  COMT  and  1MAG  are  performed  by  subroutine 
BLNNNN,  and  the  function  of  block  STOP  is  performed  by  the  main  program. 

The  executable  functional  blocks  contained  in  the  MEF  are: 

BLOCK  FUNCTION 

SOLR  Assemblage  of  distributed  load. 

L1NM  Solution  of  linear  problem  with  global  matrix  in  core. 

LIND  Solution  of  linear  problem  with  global  matrix  out  of  core. 

NLIN  Solution  of  stationary  non-linear  problem. 

TEMP  Solution  of  unsteady  problem  (linear  or  non-linear). 

VALP  Eigenvalues  and  eigenvectors. 

The  various  blocks  designed  for  execution  of  the  various  computations  have 
similar  structures  since  they  have  to: 

•  construct  element  and  load  matrices; 

•  assemble  global  matrices  and  vectors; 

•  factorize  and  solve  the  system  of  equations; 

<) 


•  output  the  results. 

Using  the  constructed  elements  and  load  matrices,  the  subroutine  element  library 
ELEMLB  is  called.  This  library  contains  subroutines  that  define  the  individual  element 
types.  The  ELEM03  subroutine,  a  thirty-two  node,  three  dimensional  isoparametic 
element  was  developed  and  added  to  the  element  library.  This  new  element  allowed  the 
solution  of  linear  elastic  structures  composed  of  homogeneous  and  isotropic  materials. 

Multiple  sample  problems  were  developed  to  fully  exercise  use  of  this  new 
element.  An  indepth  investigation  was  then  conducted  to  determine  the  limit  of 
computational  ability  using  the  newly  defined  element  to  represent  physical 
phenomenons. 


II.  DESCRIPTION  OF  THE  COMPUTER  PROGRAM 
A.  REFERENCE  ELEMENT 

To  simplify  the  analytical  expression  for  elements  of  complex  shapes  an  element 
of  reference  is  introduced.  Such  an  element  is  defined  in  an  abstract  non-dimensional 
space  with  a  very  simple  geometrical  shape.  The  geometry  of  the  reference  element  is 
then  mapped  into  the  geometry  of  the  real  element  using  geometrical  transformation 
expression. 

A  thirty-two  node,  three  dimensional  cubic  element  was  introduced  as  the 
reference  element.  The  element  has  eight  coner  nodes  and  twenty  four  mid-side  nodes 
dividing  each  edge  in  three  equal  parts  as  shown  in  Figure  2.1.  Using  this  reference 
element  we  created  the  fundamental  matrices  and  vectors  in  subroutine  called 
ELEM03,  NI03,  D03  and  B03  to  be  used  in  element  library  subroutine  ELEMLB  of 
the  MEF  program. 
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Figure  2.1  Reference  Element. 
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B.  SUBROUTINE  ELEM03 

Multi-purpose  program  MEF  is  organized  to  contain  a  library  of  one,  two,  and 
three-dimensional  elements  for  the  solution  of  problems  from  a  wide  variety  of 
disciplines.  Problems  from  the  mechanics  of  solids  and  fluids,  heat  transfer,  etc.  have 
been  solved.  For  each  element  type  nn ,  subroutine  ELEMnn  controls  the  computation 
of  all  matrices  and  vectors. 

In  element  type  03,  subroutine  ELEM03  is  organized  to  create  the  fundamental 
matrices  and  vectors  that  must  be  numerically  integrated  using  methods  of  Numerical 
Integration  described  in  and  the  computations  can  be  carried  out  by  the  following 
steps. 

1.  Operation  common  to  all  element  of  the  same  type: 

-  compute  the  weight  wf  and  the  coordinate  of  integration  points; 

-  construct  the  functions  N  (interpolation  functions),  the  function  N 
(geometrical  interpolation  functions)  and  their  derivatives 

with  respect  to  tj.  £  at  the  points  of  integration. 

2.  Operation  for  the  computation  of  matrix  [k]  of  each  element: 

-  initialize  the  matrix  [k]; 

-  for  each  point  of  integration  : 

•  construct  the  Jacobian  matrix  [J]  from  the  derivatives  with  respect 
to  x\,  (J  of  function  N  and  the  nodal  coordinates  of  the  element; 

•  construct  the  inverse  of  [J]  and  its  determinant; 

•  construct  the  derivatives  of  functions  N  with  respect  to  x,  y,  z 
starting  from  the  derivatives  with  respect  of  x\,  (J  ; 

•  construct  the  matrices  [B]  and  [D]; 

•  accumulate  into  [k]  the  values  of  [B]t[D][Bldet(J)wr  calculated  for 
each  integration  point. 

3.  Operation  required  to  compute  mass  matrix  [mj: 

-  initialize  [mj 

-  for  each  integration  point  t,T : 

•  compute  Jacobian  matrix  and  its  determinant; 

•  accumulate  the  values  of  {N}  <  N>  det(J)wr  into  matrix  [mj. 

4.  Operation  required  to  compute  consistent  load  vectors  (f): 

-  initialize  (f); 

-  for  each  integration  point  ■ 

•  compute  the  Jacobian  matrix  and  its  determinant; 


•  accumulate  the  values  of  {\'}fv  det(J)wr  into  {{]. 

5.  Operation  required  to  compute  the  residue  vecter  {r}: 

-  using  the  value  of  {1}  from  4; 

-  for  each  integration  point  £r  : 

•  compute  matrices  IB),  ID],  jj]  as  in  2  above; 

•  accumulate  the  product:  {f}  -  [B]t[D][B]{un}  wrdet(J)  into  {r}. 

6.  Operation  required  to  compute  gradients  du  at  points  of  integral: 

-  for  each  integration  point  £,r : 

•  construct  matrix  [B]  as  in  2  above; 

•  compute  and  print  gradient  {du}  =  [B){un}. 

The  subroutine  ELEM03  executes  one  operation  at  a  time  depending  on  the 
value  of  ICODE.  Control  variable  ICODE  specifies  which  element  operation  is 
desired,  and  expression  of  this  variable  as  follows: 

ICODE  =  1  initialization  of  the  characteristic  parameters  of  an  element 
(number  of  nodes,  number  of  degrees  of  liberty). 

ICODE  =*  2  operations  required  by  a  given  referance  element  which  are 
independent  of  the  real  geometry;  construction  of 
interpolation  function  N  and  their  derivative  with 
respect  to  4  at  the  points  of  integration. 

ICODE  *  3  construction  of  matrix  Ik]  in  array  VKE. 

ICODE  =  4  construction  of  tangent  matrix  needed  for  non-linear 
problems  in  array  VKE. 

ICODE  =  5  construction  of  mass  matrix  [m]  in  array  VKE. 

ICODE  =  6  computation  of  residual  vector  {r}  in  array  VFE. 

ICODE  =  7  computation  load  vector  {f}  in  array  VFE. 

ICODE  =  8  computation  and  printing  of  gradients  {<3u}. 


C.  CODING. 

1.  Evaluate  coordinates,  weights,  functions  N  and  their  derivatives 

Integration  formular  for  numerical  integration  is  the  following  form. 


where 


(eqn  2.1) 


-  are  the  coordinates  of  integration  point  T  in  £,,  H.  £  system 
coordinate  corresponding  to  weight  Wj  ; 


13 


-  wj  are  the  weights  corresponding  to  integration  point  number; 

-  IPG  are  the  total  number  of  integration  points. 

A  choice  of  2,  3  or  4  integration  points  by  dimension  can  be  made  in 
subroutine  GAUSS  [Ref.  l:p.  265],  giving  respectively  8,  27  and  64  integration  points, 
the  weights  corresponding  to  the  integration  points  and  their  coordinates.  They  are  in 
the  array  called  IPG,  VCPG  and  VKPG  respectively. 

Subroutine  NI03  create  the  array  VNI  that  contains  the  shape  function  Nj 
and  their  derivatives  with  respect  to  £  system  coordinate  )  as  Figure  2.2. 


Figure  2.2  Block  Diagram  for  the  Shape  Function  and  their  Derivative. 

2.  Computation  of  stiffness  matrix  [kj  of  each  element. 

The  explicit  equation  of  element  stillness  matrix  is  following: 

lit 

[k]  =  J  J  J  [B|t[Dj[B!det<J)d5  dn  (eqn  2.2) 

-i  -i  -l 

where 

-  [B[  is  the  linear  strain  matrix; 

-  [B)c  is  the  transpose  matrix  of[B); 

-  [D]  is  the  matrix  of  elastic  constants  for  an  isotropic  material. 

The  stillness  matrix  that  has  the  block  diagram  as  Figure  2.3  is  the  integration 
of  product  of  the  linear  strain  transpose  matrix,  the  element  property  matrix  and  the 
linear  strain  matrix  over  the  volume  of  the  reference  element. 


[k°]  =  (B]t[DJ[B]det(J) 


(eqn  2.4) 


Subroutine  JACOB  (Ref.  l:p.  63],  compute  the  Jacobian  matrix,  its  inverse 
and  its  determinant.  The  Jacobian  matrix  as  Figure  2.4  is  obtained  as  the  product  of 
two  matrices,  one  containing  the  derivatives  of  the  geometrical  transformation 
functions  with  respect  to  the  space  of  the  reference  element,  and  the  other  containing 
the  real  coordinates  of  the  geometrical  nodes  of  the  element. 

<  xyz  >  =  <N(£)>[{xn}  {yn}  {zn}]  (eqn  2.5) 

{xn},  {yn},  {zn}  being  the  geometrical  node  coordinates.  The  Jacobian  matrix  is  figure 
2.4 


Figure  2.4  The  Jacobian  Matrix. 

The  inverse  of  Jacobian  matrix  [J]**  as  Figure  2.5  and  its  determinant  are 
computed  by  the  method  discribed  on  (Ref.  l:p.  44], 

Subroutine  DNIDX  [Ref.  l:p.  64],  computes  the  derivatives  of  the  shape 
function  i\-  with  respect  to  the  coordinate  system  of  the  real  element  using  the  product 
on  Figure  2.6  . 

Subroutine  B03  creates  the  matrix  as  Figure  2.7  that  contains  the  strain 
components  at  all  direction  of  each  node  in  the  element  using  output  components  of 
subroutine  DNIDX  and  rearrange  the  new  array  called  VBE. 

Subroutine  D03  compute  the  stress-strain  matrix  (VDE)  as  Figure  2.S 
(Ref.  2:p.  110]  for  isotropic  materials  (  E  =  Young's  Modulus,  v  =  Pisson's  Ratio  ). 
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Figure  2.5  The  Inverse  of  Jacobian  Matrix. 


Figure  2.6  The  derivative  of  shape  function  w.r.t.  x,  y,  z. 

Subroutine  BTDB  construct  the  stiffness  matrix  by  adding  the  product  of  the 
transpose  matrix  VBE,  the  stress-strain  matrix  VDE  and  the  matrix  VBE  of  every 
integration  points. 

[k]  =  (B]t[DJ[B]widet(J)  (eqn  2.6) 


3.  Computation  The  Mass  Matrix  [m). 

The  element  mass  matrix  has  the  block  diagram  as  Figure  2.9  and  the  explicit 
equation  as  following: 


[m]  -  Jj’j  jW[N]det(J)  d$  dq 


(eqn  2.') 


and  numerical  integration  form  is: 


•  •  "l,«  *  # 


0  _  0  0  H, 


H.  *, 

ii>  *.r 


«  o 


•  "».r  * 


"ji.r  "m,«  0 
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J*.»  M.r 


'j*t«  #  "u,« 


Figure  2.7  The  Array  VBE. 
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Figure  2.8  The  Array  VDE. 

[m]  =1?wi[N]t[Nldet(J) 
i=l 


(eqn  2.8) 


-  [N]  is  the  shape  function  matrix  Figure  2.10  that  was  created  by 
subroutine  NI03; 

-  [N)t  is  the  transpose  matrix  of  [N]; 

-  w-  are  the  weighting  coefficients  corresponding  to  each  integration 
point; 

Obviously  look  at  the  product  of  matrices  [N]land[N]  as  Figure  2.11  is  a 
symmetric  matrix.  3y  convention,  the  mass  matrix  [mj  contains  upper  half  components 
and  diagonal  components  of  this  matrix. 


Figure  2.11  The  Product  of  Matrices  [N]1  and  [N]. 


fyx*  fyy,  fvz  are  the  force  per  unit  volume  in  direction  x,  y,  z. 

5.  Computation  the  residue  array. 

The  element  residual  array  has  the  block  diagram  as  Figure  2.13  and  defined 
by: 


{r}  =  (0  -  [k]{un} 


(eqn  2.1 1) 


where 


-  {r}  is  the  residue  vector; 

-  {f}  is  the  element  load  vector; 

-  [k]  is  the  element  stiffness  matrix; 

-  {un}  is  the  nodal  values  vector. 

6.  Computation  of  the  gradients  and  stresses  at  the  points  of  integration. 
For  each  integration  point: 

-  compute  strain  as  Figure  2.14.  [Ref.  3:p.  31] 

or  compacted  form 


{c}  =  [B]{Ui} 


(eqn  2.12) 


fi 


,-^j  _  -  compute  stress  as  Figure  2.15.  [Ref.  3:p.  30] 

or  compacted  form 


.•  _v_V  . 


-*»  Ji  A  Ja  A  . 
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Figure  2.14  The  Strain-Node  Value  Relation. 
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Figure  2.15  The  Stress-Strain  Relation. 
{<Tj}  =  [D]^} 


(eqn  2.13) 


compute  coordinate  of  each  integration  points  on  the  real  element  as 
Figure  2.16. 
or  compacted  form 


{^}  -  [N]{xni} 


(eqn  2.14) 


Using  xn,  yn,  zn  .where  N  =  N(£j,  r[v  ^)  and  ^  ^  are  the  gauss  points, 

the  coordinate  of  thirty  two  nodes  on  the  real  element  that  existed  in  the  array 
VCORE  .  All  the  integration  points  are  evaluated  in  the  reference  element.  Using  the 
product  of  the  transfer  matrix  ( the  shape  function  matrix  N)  and  the  array  VCORE,  we 
can  evaluate  the  coordinate  of  the  integration  points  on  the  real  element. 


Figure  2.16  Evaluation  of  the  Coordinate  of  the  Integration  Points  on  the  Real  Element. 


Note  that  subroutine  ELEM03  executes  one  operation  at  a  time  depending  on 
the  value  of  ICODE.  For  preserving  the  memory  locations,  some  of  the  operations 
have  been  performed  using  the  same  array  name  as  VKE  in  both  the  stiffness  matrix 
[k]  and  the  mass  matrix  [mj.  The  same  thing  has  been  done  with  array  VFE  used  for 
the  residual  vector  {r}  and  the  load  vector  {!}. 


III.  THE  SOLUTION  OF  A  SIMPLE  PROBLEM  AND  DISCUSSION  OF 

RESULT 


In  this  .chapter  the  preparation  of  input  data  for  the  computer  program  are 
described.  A  simple  problem  was  developed  and  the  investigation  was  conducted  to 
determine  the  ability  of  this  cubic  element. 

A.  ENTRY  AND  EXECUTION  FUNCTIONAL  BLOCKS 

MEF  has  specialized  functional  blocks  for  the  entry,  verification  and 
organization  of  the  data  required  to  define  a  problem.  Block  COOR  reads  the  nodal 
coordinates  and  the  number  of  degrees  of  freedom  of  each  node,  it  also  provides 
automatic  node  generation.  Block  COND  reads  the  boundary  condition.  Block  PREL 
reads  element  properties  if  required  for  element  type  being  used.  Block  SOLC  reads 
the  concentrated  loads.  Block  ELEM  reads  the  element  connectivities;  it  also  reads 
element  group  information  when  more  than  one  element  type  is  used,  if  elements  have 
different  properties.  This  block  provides  automatic  element  generation. 

Other  function  blocks  of  MEF  for  the  execution  of  particular  finite  element 
computations  use  the  data  base  constructed  by  entry  blocks  and  augment  it  by  their 
results.  Block  LINM  assembles  and  solves  a  linear  system  of  equations  residing 
in-core.  Block  LIND  is  similar  to  the  block  LINM  but  the  system  of  equations  resides 
out-of-core  in  a  mass  storage  device.  And  block  STOP  terminates  execution  of  the 
problem. 

MEF  provides  various  levels  of  output.  The  quantity  of  output  desired  from  a 
given  block  is  controlled  by  a  parameter  on  the  block  calling  card,  described  in  detial 
[Ref.  l:pp.  440-447],  which  ranges  from  0  (the  assumed  value)  to  4.  The  default  value 
provides  all  the  information  needed  to  verify  the  input  stream  and  obtain  the  desired 
answers  while  the  value  1  thru  4  provide  various  level  of  verbosity. 

B.  STRUCTURE  MODELING 

The  first  step  in  applying  a  finite  element  solution  to  the  problem  is  the  selection 
of  an  appropriate  mesh.  In  many  case  an  extrapolation  of  the  results  will  be  required 
and  hence  more  than  one  mesh  will  have  to  be  selected.  The  mesh  size  solution  does 
not  follow  any  predetermined  rules  and  will,  in  general,  depend  on  the  nature  of  the 
problem  and  the  judgement  of  the  analyst.  An  abitrarv  numbering  convention  is 
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adopted  for  the  nodal  points  of  the  entire  structure  (in  distinction  with  the  numbering 
convention  for  the  element  nodes  shown  in  Figure  3.1.  A  second  numbering  scheme  is 
also  needed  for  the  elements  of  a  mesh. 

In  order  to  show  the  function  of  the  program,  as  well  as  some  observations 
affecting  the  use  of  it,  the  solution  of  a  simple  problem  is  presented  in  this  section. 
The  problem  selected  is  that  usually  presented  in  classical  texts  of  strength  of  materials 
as  a  cantilever  beam  of  uniform  cross  section  subjected  to  a  concentrated  load  at  the 
end  of  the  beam.  The  model  used  for  this  problem  is  shown  in  Figure  3.1.  Nodes  1,  2, 
3,  4,  5,  6,  7,  8,  9,  10,  11  and  12  (here  arbitrarily  numbered)  are  constrained  to  zero 
displacement  in  all  directions.  The  concentrated  load  at  the  end  of  the  cantilever  beam 
is  considered  as  the  consistent  load  and  are  shown  in  Figure  3.2. 

Four  different  meshes  which  each  mesh  has  the  thickness  9,  0.9,  0.09,  0.009 
inches  and  has  the  concentrated  loads  1200,  1.2,  1.2e-3  and  1.2e-6  pounds  were 
employed  in  order  to  show  the  variation  of  results  with  mesh  size.  Numerical  results 
for  the  maximum  displacement  at  the  tip  of  the  beam  are  shown  in  Tables  I  and  II. 


Figure  3.2  The  consistent  load  at  the  end  of  the  beam. 

C.  DISCUSSION  OF  RESULTS 

Elementary  beam  theory  gives  a  displacement  of  0.0631  inches,  however  this 
beam  theory  neglects  shear  deflection  and  three  dimensional  elasticity  does  not  do  so. 
Observation  of  the  results  presented  here  reveal  that  mesh  refinement  lead  to  improved 
results  which  eventually  will  converge  to  a  certain  value.  The  thirty  two  nodal  points 
brick  is  subject  to  numerical  bad  conditionary  when  used  with  an  adverse  slenderness 
ratio.  This  ratio  being  defined  as  the  ratio  of  the  maximum  element  dimension  over 
the  minimum  dimension.  For  example  in  the  numerical  examples  treated  the 
slenderness  ratio  for  the  one  and  eight  elements  representation  of  the  four  beams 
analysis  are  shown  in  the  Tables  I  and  II.  Results  are  acceptable  for  the  first  three 
beams  in  each  case.  However  a  computation  of  approximate  values  of  zero,  first 
invariant  of  the  element  stiffness  as  well  as  condition  number  are  produced  to 
investigate  the  results. 

As  the  slenderness  ratio  increases  the  numerical  of  the  conditioning  of  the 
stiffness  matrices  become  so  bad  that  no  significant  digits  can  be  expected  out  of  the 
solution  for  displacements.  The  same  conclusion  is  arrived  with  a  mesh  of  eight 
elements.  In  such  a  case  the  slenderness  ratio  varies  from  1.6  to  1666.6.  For  the 
thickness  of  0.09  the  slenderness  ratio  was  adequate.  To  reduce  the  ratio,  more 
element  must  be  used  and  tests  with  20  and  40  elements  were  conducted.  The  twenty 
elements  mesh  gives  of  0.6  to  666.6.  In  the  use  of  a  mesh  with  fourty  elements  the 
largest  dimension  become  6.0  and  must  therefore  be  reduced  to  3.0  as  3.0  is  the  value 


of  120/40.  Runs  with  such  a  mesh  were  then  performed  to  confirm  our  results.  On  the 
basis  of  all  these  runs  and  the  value  of  the  condition  number  of  these  matrices  it  seems 
that  a  slenderness  ratio  of  approximately  150  could  still  be  used  to  give  results  with  6 
significant  digits  in  the  answers,  however  further  examples  must  be  treated  before  a 
conclusion  could  be  reached. 

These  are  observations  essentially  made  on  the  results  of  a  simple  problem.  Thus 
no  firm  rules  regarding  the  use  of  the  newly  element  can  be  established  and  further 
experimention  with  different  problems  has  to  be  conducted  for  this  purpose. 


TABLE  I 

THE  DISPLACEMENT  FOR  THE  ONE  ELEMENT 


Sloodornoss 

Ratio  . 

mm 

Thickness 

Zaro 

13 

2>, 

^msx 

^min 

Nuabor 

Displacement 

9 

13.3 

- 

0.24«-3 

-0. 5*0-3 

7.27all 

.  3. 72elO 

3.l*a** 

I.I80O6 

-0.0381 

0.9 
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O.76o-S 

-0. 16a-* 
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0.09 
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TABLE  II 

THE  DISPLACEMENT  FOR  THE  EIGHT  ELEMENTS 


Thickness 
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Ratio 
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APPENDIX 

ELEM03  SUBPROGRAM  LISTING 


The  listing  of  subprogram  ELEM03  is  provided  below.  It  includes  four 
subroutines:  NI03,  D03,  B03  and  BTDB  respectively. 


It  It  it  It  it  it  it 


SUBROUTINE  ELEM03( VCORE, VPRNE , VP  REE , VDLE, VKE, VFE) 

***********************k*x*********x****x********' 


it  it  it  it  irk  it  it  it  it  it  Hit  it  it  it 


32  NODES  HEXAHEDRON  ELEMENT  FOR  3  DIMENSIONAL  ELASTICITY 
EVALUATE  ELEMENT  INFORMATIONS  ACCORDING  TO  ICODE  VALUE 
IC0DE=1  ELEMENT  PARAMETERS 

IC0DE=2  INTERPOLATION  FUNCTIONS  AND  GAUSS  COEFFICIENTS 
IC0DE=3  STIFFNESS  MATRIX 
IC0DE=4  • 

IC0DE=5  MASS  MATRIX 
IC0DE=6  RESIDUALS 
IC0DE=7  SECOND  NUMBER 
IC0DE=8  EVALUATE  AND  PRINT  STRESSES 
ELEMENT  PROPERTIES 

VPREE(l)  YOUNG'S  MODULUS 
VPREE(2)  POISSON'S  COEFFICIENT 
VPREEC3  i  SPECIFIC  MASS 

********************************************************************** 

IMPLICIT  REAL*8(A-H ,0-Z) 

COMMON  /COOR/NDIM 
COMMON  /ASSE/NSYM 

COMMON  /RGDT/IEL, ITPE. ITPE1.IGRE, IDLE, ICE, IPRNE, IPREE, INEL, IDEG, 

1  IPG  ICODe, IDLES, INEf.0, IP60 

COMMON  /ES/MAMR,tiP 

DIMENSION  VC6RE( 1 ) , VPRNE( 1 ) , VPREE( 1 ) , VDLE( I ) , VKE( 1 ) , VFE( 1 ) 


CHARACTERISTIC  DIMENSIONS  OF  THE  ELEMENT 


DIMENSION  VCPG( IPG) ,VKPG(NDIM*IPG) .VDE1( IMATD**2) 

DIMENSION  VCPG(270KPG($1).VDE1(36) 

DIMENSION  VBEuMATO*IDLE).v6E(IMATD**2),VJ(NDIM**2),VJ1(NDIM**2) 
DIMENSION  VBEC5761, VDEC365 ,VJ  9),VJ1(9) 

DIMENSION  VNIX(INEL*NDIM0NI(C  1+NDIM)*INEL*IPG) , IPGKED(NDIM) 
DIMENSION  VNIX(96),VNI(34$6),IPGKED(3) 

—  DIMENSION  OF  MATRIX  D,  NUMBER  OF  G.  P. 

DATA  IMATD/6/. IPGKED/3 , 3. 3/ 

DATA  ZERO/O. 0£)0/, DEUX/2.  60/,X05/0. 5D0/,RADN/.  572957795130823D2/ 
DATA  EPS/1. D-6/ 

DATA  NNNNNN/0/ 


CHOOSE  FUNCTION  TO  BE  EXECUTED 


GOTO  (100, 200, 300, 400 ,500 ,600 ,700 ,800), ICODE 
100  IDLE0=96 
INEL0=32 
I PG0=2  7 
RETURN 

*********************************************************************** 

EVALUATE  COORDINATES .WEIGHTS , FUNCTIONS  N  AND  THEIR 
DERIVATIVES  AT  G. P. 

*********************************************************************** 


200  CALL  GAUSS( IPGKED ,NDIM, VKPG ,VCPG, IPG) 
IFfM.LT.2)  GOTO  2^0 


.vs-.-;  ;vw 
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2000  FORMA^C/ii?^°^AUSS  POINTS 1 /10X , 1 VCPG * ,25X, 1 VKPG 1 ) 
10=1 

DO  210  IG=1. IPG 

Il=IO+NDIM-i  ,  ,  v 

210  io-io1nd1m010)  VCPQCIGJ.CVKPGCD.I.IO.Il) 

2010  FORMATflX, F13.  9.5X.3F13.  9) 

220  CALL. NI03( VKPG, #NT< 

IF(M.LT.  2)  retOrn 


I1=4*INEL*IPG. 
u( 


2020  FORMAT(I/1X^^F^nSt¥o^IAkID— DERIVATIVES 1  /( 1X,8E12.  5)) 

RETURN  " 

^*******iiMit*******)inif*****jlf********#i*r***#***>l(*lk<rjlr***ilf*jlt*jifjinir*#jlf*#HfHr>lfilt*>lr*>if* 

C  EVALUATE  ELEMENT  STIFFNESS  MATRIX 

0********************************************************************** 

C - INITIALIZE  VKE 

300  NINI=4656 

IF(NSYM.NE.O)  NINI=9216 
DO  310  1=1 .NINI 
310  VKE(I)=ZERd 

C - --  FORM  MATRIX  D 

CALL  D03( VPREE.VDE) 

C - LOOP  OVEft  THE  G.  P. 

I1=1+INEL 
DO  330  IG=1.IPG 

C -  EVALUATE  tHE  JACOBIAN. ITS  INVERSE  AND  ITS  DETERMINANT 

c _ CALLJACOB^VNiyi^pftE,NDIM,INEL,VJ,VJl,DETJ) 

C=VCPG(IG)*DETJ 
DO  320  1=1.36 
320  VDEl(I)=VD£m*C 
: - -  PERFORM  MATRIX  B 


NDIM, INEL, VNIX) 
IDLE, IMATD.NSYM) 


CALL  DNIDX  , 

CALL  B03CVN] 

CALL  BTDB  ( 

330  I1=I1+4*INE 

RETURN 
400  CONTINUE 
RETURN 

Q*  *  *  *  it  *  *  *  **  it  *  *  ★  *  *  *  ★  it  ft  *  ★  *  *  ★  *  **  *  *  ★  *  ft  ft  ft  ft  ★  ft  ft  *  *  ft  Hr  *  ★  *  *  ft  it  it  *  **  it  it  *  *r  ★*  *  *  ★ 

C  EVALUATE  THE  MASS  MATRIX 

Q********************************************************************** 

500  NINI=4656 

IF(NSYM.NE.O)  NINI=9216 
DO  510  1=1 .NINI 
510  VKE(I)=ZERO 

C  LOOP  OVER  THE  G. P. 


C 

C- 

C 


I  DIM 1=ND IM— 1 
I  DEC  L=(ND IM+1 ) * I NEL 
I1=1+INEL 
12=0 

DO  550  IG=1 , IPG 

CALL  JACOB(VNI( II), VCORE, NDIM, INEL, VJ,VJ1,DETJ) 
D=VCPG( IG)*DETJ*VPREE( LL) 


ACCUMULATE  MASS  TERMS 


IDL=0 

DO  540  J=1 , INEL 
JJ=I2+J 

J0=1+IDL*( IDL+l)/2 
DO  530  1=1, J 
11=12+1 

C=VNI( II)*VNI(JJ)*0 
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VKE(J0)=VKE(J0)+C 
Jl=J0+lDL+2 
DO  520  I 1=1 . IDIM1 
VKE(J1)=VKE(J1)+C 
520  Jl=Jl+lDL+3 
530  J0=JQ+NDIM 
540  I DL=I DL+NDIM 
1 1=1 1+IDECL 
550  I2-I2+IDECL 
RETURN 

C  EVALUATE  THE  ELEMENT  RESIDUAL 

Q****************************************************************** 

C - FORM  MATRIX  D 

600  CALL  D03(VPREE,VDE) 

C -  INITIALIZE  THE  RESIDUAL  VECTOR 

C 

DO  610  I D=1 . IDLE 
610  VFE(ID)=ZER6 

- L00p  OVER  THE  G.  P. 

I1=1+INEL 
DO  640  IG=1 .IPG 

-  EVALUATE  THE  JACOBIAN 

CALL  JAC0B( VNI( II) . VC0RE.NDIM. INEL. VJ , VJ1 ,DETJ) 

- EVALUATE  FUNCTIONS  D(Nl)/D(X) 

CALL  DNIDX(VNI( II) , VJ1 ,NDIM, INEL, VNIX) 

-  EVALUATE  STRAINS  AND  STRESSES 

EPSX=ZERO 
EPSY=ZERO 
EPSZ=ZER0 
GAMXY=ZERO 
GAMYZ=ZERO 
GAMZX=ZERO 
ID=1 

DO  620  IN=1 , INEL 
UN=VDLE( ID) 

VN=VDLE(ID+1) 

WN=VDLEfID+2; 

Cl=VNIXf IN) 

IN1=IN+INEL 
C2=VNIX( INI) 

IN2=IN1+INEL 
C3=VNIX( IN2) 

EPSX=EPSX+Cl*UN 
EPSY=EPSY+C2*VN 
EPSZ=EPSZ+C3*WN 
GAMXY=GAMXY+C1*VN+C2*UN 
GAMYZ=GAMYZ+C2*WN+C3*VN 
GAMZX=GAMZX+C1*WN+C3*UN 
620  ID=ID+3 

C1=VCPG( IG)*DETJ 
C2=VDE(2)*C1 
C3=VDE(22)*C1 
C1=VDE(1)*C1 

SIGX=Cl*tPSX+C2*EPSY+C2*EPSZ 
SIGY=C2*EPSX+C1*EPSY+C2*EPSZ 
SIGZ=C2*EP$X+C2*EP$Y+C1*EPSZ 
TAUXY=C3*GAMXY 
TAUYZ=C3*GAMYZ 
TAUZX=C3*GAMZX 

. —  FORM  THE  RESIDUAL 


10=1 

DO  630  IN=] 
C1=VNIX(IN 
IN1=IN+INEL 
C2=VNIX( INI' 
IN2=IN1+INE( 


INEL 


VFE 


630 

640 


VFE 

yfE 


VNIX 

IN2) 

(ID)= 

VFEflC 

n+ci* 

i  ID+1 

)=VFE( 

ID+1) 

flD+2 

)=VFE| 

ID+2) 

ID+3 

+C2*SIGY+C1*TAUXY+C3*TAUYZ 

+C3*SIGZ+C2*TAUYZ+C1*TAUZX 


I1=I1+4*INEL 
RETURN 

£********************************************************************** 
C  EVALUATE  BODY  FORCES, FX.FY.FZ  PER  UNIT  VOLUME 

(FOR  GRAVITY  FX=0 ! FY^O , ^Z=- VPREE( 3 ) 
************************************ 


£************* 

700  FX=ZERO 
FY=ZERO 
LL=3 

FZ=-VPREE(LL) 
‘  710  1=1.96 


*********************** 


710 


‘*INEL 

rj 


DO  710  1=1  « 

vfeJi)=zer6 

I DECL=( NDIM+1 ) ” _ 

DO  730  IG=1.IPG 

CALL  JACOBIAN  I Fll+INEL) , VCORE.NDIM,  INEL,  VJ ,  VJ1  ,DETJ) 
DX=yCPG(IG)*DETJ  '  ' 

DY=DX*FY 
DZ=DX*FZ 
DX=DX*FX 
12=11 
13=1 

DO  720  IN=1, INEL 
VFE( I3)=VFE( 13 j+DX*VNI(I2) 
*--i-/l)=VPE(l3+l)+DYvVNf(I2 


720 

730 


VFE( 13+1 
VFE  ■’  ‘ 
12= 

13= 

11= 


B) 


[I3+2)=VFE(I3+2)+DZ*VNl 
[3+3 

Il+IDECL 
RETURN 

Q******************************************************************* 

C  EVALUATE  AND  PRINT  STRESS  AT  G.  P. 

Q******************************************************************* 

800  WRITE(MP,2080)  IEL 

2080  FORMAT(/;Y  STRESSES  IN  ELEMENT’  1 5/ 1  ,  P. G.  ,6X, 1 X’  ,  11X , ’ Y1 , 11X, 
1  Z’ ,9X , 1 SIGX 1 ,8X,  SIGY  ,8X,  SIG2’ ,/X,  TAUXY  7X * ’ TAOYZ ’  7X, 

2 '  TAClZX '/) 

c 

C - FORM  THE  MATRIX  D 

C 

CALL  D03( VPREE, VDE) 

C - LOOP  OVER  THE  G.  P. 

I 1=1+INEL 
12=0 

DO  820  IG=1 . IPG 

C— . EVALUATE  THE  JACOBIAN 

CALL  JACOBI VNI(I1). VCORE.NDIM. INEL, VJ.VJl.DETJ) 

C - EVALUATE  FUNCTIONS  D(NI)/£KX) 

CALL  DN IDX( VN 1(11), VJl.NuIM, INEL, VN IX) 

C . COMPUTE  STRAINS  AND  COORDINATE  AT  G.  P. 

C 

EP$X=ZERO 

EPSY=ZERO 

EPSZ=ZERO 

GAMXY=ZERO 
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(N=1,INEL 


810 


GAMYZ=ZER0 
GAMZX=ZERO 
X=ZER0 
Y=ZER0 
Z=ZER0 
ID=1 
DO  810 
UN=VDLE( 

VN=VQLE' 

WN=VDLE 
XN=VC0RL 
YN=VC0RE 
ZN=VC0RF 
C1=VNIX( _  , 

IN1=IN+:  NEC 
C2=VNIX( INI) 
IN2=IN1+INEL 
C3=VNIX(IN2) 

IN1=IN+I2 

C4=VNI( INI) 

EPSX=EPSX+Ul*UN 

EPSY=EPSY+C2*VN 

EPSZ=£PSZ+C3*WN 

GAMXY=GAMXY+C1*VN+C2*UN 

GAMYZ=GAMYZ+C2*WN+C3*VN 

GAMZX=GAMZX+C1*WN+C3*UN 

X=X+C4*XN 

Y=Y+C4*YN 

Z=Z+C4*ZN 

ID=ID+3 


COMPUTE  THE  STRESSES 


SIGX=VDE(1)*EPSX+VDE(2)*EPSY+VDE(2)*EP$Z 
SIGY=VDE(  2 j*EPSX+VDE(  ll*EPSY+VDEf 2 j*EPSZ 
SIGZ=VDE(2)*EPSX+VDE(2)*EPSY+VDE( 1;*EPSZ 
TAUXY=VDE(22)*GAMXY 
TAUYZ=VDE(22)*GAMYZ 
TAUZX=VDE(22)*GAMZX 

WRITEfMP ,2090)  IG,X.Y,Z,SIGX,SIGY,SIGZ,TAUXY,TAUYZ,TAUZX 
2090  FORMATflft, I5.9E12.5j 
I2=I2+4*II\IEL 
820  1 1=1 1+4*INEL 

RETURN 
END 

SUBROUTINE  NI03(VKPG, VNI) 

*************************k***k***************************************** 

TO  EVALUATE  THE  SHAPE  FUNCTIONS  N  AND  THEIR  DERIVATIVES  W. R. T. 
KSI.ETA.DZETA 
INPUT 

VKPG(NN)=  COORDINATES  OF  POINTS  IN  KSI.ETA.DZETA 
OUTPUT 

VNI  =  THE  SHAPE  FUNCTION  N  AND 

THE  DERIVATIVE  OF  SHAPE  FUNCTION  W.  R.  T.  KSI.ETA.DZETA 
*********************************************************************** 

IMPLICIT  REAL*8  (A-H.O-Z) 

DATAN«M/^jb4°3):S=i;3liVKPG(81),VN(32'3)'RN(32) 

*  /-l.  £)0:-l.  Dd.-i.  DO  -. 533333333333333,-1.  DO, -1.  DO, 0.  3333 

*33333333333,-1.  0(5  .-1.  0(5,1.  DO  -1.D0.-1.  DO ,  1.  DO . 3333^3333353333 A 
*-l. DO , 1. DO,. 333335333335333,-1. 00,1.00,1. DO, -1. DO,.  333333333333533 
*1.  DO  -1.  0(5,-.  333333333333353,1.  D(5,-l.  (50,-l.  DO ,  1.  (50,-1.  D0,-1.  DO, 
*(5.  333533333533333,-1.  DO ,-l.  DO  -.  335333335333335,-1.  £>0,-1.  £)0  -1.  Do, 
*-. 333333333333333,1. 00,-1. D0,:. 333333333333333,1. D0,1. DO,-. 5333335 
*33333333,-1.  00,1.  £)0,-.  533333533333333,-1.  D0,-1.  DO,.  533335333333333 
DO,:L  D0,.j53333533333333,l.  D0,1.  0(5, ., 333533333533333 ,-l.  00,1.  CO 
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*1.  DO, 1.  DO.-.  333333333333333.1. DO. 1.  DO.-l.  DO.l. DO, 1.  DO,-l.  DO,.  3333 
*33335333353, l.DO.-l. DO,-. 335333353333533,1.00/ 


DO  1  1=1,32 

WRITE?*, $9)  I  ,(CORRF( I ,J) .J=l,3) 
FORMAT  (2X, I2’3(2X,E12.  6)) 

DO  10  NN=1 ,81,3 
X=VKPG(NN) 

Y=VKPGfNN+l) 

Z=VKPG(NN+2l 
PRINT*'x.Y.Z 
-  At  dote 

-  NODE  #  :  1,4,7,10,21,24,27,30 


DO  100  1=1.2 

II=20*(I-1)+1 

IT=II+9 

DO  100  J=I I . IT , 3 
Xl=CORRF(J,i) 

Y1=C0RRF(J  2) 

Zl=CORRFfJ  3) 

CF1=9.D0>64.D0 
CF2=19.  DO/9.  DO 

RN(J)=CF1*(1. D0+X*X1)*(1. D0+Y*Y1)*( 1.  D0+Z*Z1)*(-CF2+X*X+Y*Y+Z*Z) 
VN(J,l)=CFl*f 1.  DO+Y*t 1)*( 1.  D0+Z*Zl)*(Xl*(-CF2+3. DO*X*X+Y*Y+Z*Z) 

*  +2.D0*X)  v 

VN(J ,2)=CF1*( 1. D0+X*X1)*( 1. D0+Z*Zl)*(Yl*(-CF2+X*X+3. DO*Y*Y+Z*Z) 

*  +2. DO*Y) 

VN(J ,3)=CF1*^1.  D0+X*X1)*( 1.  D0+Y*Yl)*(Zl*(-CF2+X*X+Y*Y+3.  DO*Z*Z) 


AT  MIDSIDE 

NODE  #  :  2,3,8,9,22,23,28,29 


DO  200  K=l, 
IK=20*(K-1) 
DO  200  1=1, 
II=6*(I-1)+ 
IT=I 1+1 


=12 

l)+2+-IK 


DO  200  J=II ,  IT 
X1=C0RRF( J , 1 ) 
Y1=C0RRF(J  2' 
Z1=C0RRF(J ’ 3) 
CF1=81.  Du/64.  DO 
CF2=1. DO/9. DO 
RN(J)=CF1*(1. DO 
VNfJ.l^CFWl. 


VN(J  2)=CF1*Y1 
200  VN(j;3)=CFl*Zl 


1,  DQ+Y*Y1 


D0-X*X)*(CF2+X*X1)*( 1.  D0+Y*Y1)*( 1.  [ 
Cl. DO+Y*Yl)*(l. D0+Z*Z1 )*(Xl-2.  DO^X*C 
Yl*(l.  D0-X*XWCF2+X*X1)*(1. DO+Z*Zi; 
Zl*(l.  D0-X*X)*(CF2+X*X1)*(1.D0+Y*Y1 


1.  D0+Z*Z1) 

X*CF2-3.  D0*X*X*X1) 


AT  MIDSIDE 

NODE  #  :  5,6,11,12,25,26,31,32 


DO  300  K=1 


IK=20*(K-l) 

DO  300V 1=1,2 
I I=6*(I-l)+5+IK 
IT=II+1 
DO  300  J=II. 

X  1=C0RRF( J , 1 
Y1=C0RRF(J  2 
Z1=C0RRF(J  3 
RN(J)=CFl*{ 1 
VN(J , l)=CFl* 


)=CFl*{ 1. 00+X*Xl)*(l.  DO-Y*Y 
,  1  )=CFi*Xl*( 1.  DO-Y^Y)*( CF2+ 


1. D0+Z*Z1 


)0+Z*Zl) 


VNfJ  2)=CF1*(1. D0+X*X1)*  1. D0+Z*Z1)*(Y1-2.D0*Y*CF2-3.D0*Y*Y*Y1) 
300  VN(j;3)=CFl*Zl*(l.D0+X*Xl)*(l.D0-Y*Y)*(CF2+Y*Yl) 


AT  MIDSIDE 

NODE  #  :  13,14,15,16,17,18,19,20 


nononoono 


DO  400  J=13.20 
X1=C0RRF( J , 1 
Y1=C0RRF£J  2 
Z1=C0RRF( J  3 


RN(J)=CFWl.D0+X*Xl)*(l.p0+Y*Yl)*(l.D0-Z*Z)*(CF2+Z*Zl) 

VNrj,l)=CFl*Xl*(l.D04Y*YlWl.D0-Z*zWCF2+Z*Zl) 


VNfJ  2)=CFl*Yl*f 1.  D0+X*Xl)*fl.  D0-Z*ZW CF2+Z*Zl! 

400  VN(J!3)=CFl*fl. D0+X*X1)*(1.  DQ+Y*Yl)*(Zl-2. D0*Z*CF2-3. D0*Z*Z*Z1) 


DO  4l0  L=l,3i 

jj=jj+1=^  L 


410 

DO  420  K=1 ,3 
DO  430  1=1  32 
VNI(JJ)=VN( I ,K) 

430  JJ=JJ+1 
420  CONTINUE 
10  CONTINUE 
RETURN 
END 

SUBROUTINE  D03(VPREE, VDE) 
****************************** 


***************************************** 


TO  FORM  MATRIX  D  (  3  DIMESIONAL  ELASTICITY) 

INPUT 

VPREE  =  ELEMENT  PROPERTY 

VPREE(l)  =  YOUNG'S  MODULUS 
VPREE(2)  =  POISSON'S  RATIO 

OUTPUT 

VDE  =  MATRIX  D 

*********************************************************************** 


IMPLICIT  REAL*8 
DIMENSION  VPREEf 
DATA  ZER0/0. 000/ 
E=VPREE(1) 

LL=2 

X=VPREE(LL) 
Cl=E*(UN-X)/(( 
C2=C1*X^UN-X) 


,Ufsl/l.  OD07, 


DEUX/2. 0D0/ 


UN+X)*(UN-DEUX*X)) 


C3=C1*(0N-DEUX*X)/(DEUX*(UN-X)) 
DO  10  J=1 , 36 
10  VDE(J)=ZEftO 


VDEf 1)=C1 
VDEf2)=C2 
VDEf3)=C2 
VDEf 7)=C2 
VDE( 8 j=C 1 
VDEf 9}=C2 
VDEf 1 3 )=C2 


VDEf 14)=C2 
VDEf 15  )=C1 
VDE(22)=C3 
VDE( 29)=C3 
VDE(36)=C3 
RETURN 
END 

SUBROUTINE  B03f VNIX  INEL, VBE) 
****************************** 


***************************************** 


Q************* 

C  TO  FORM  MATRIX  B  (3  DIMEMSIONAL  ELASTICITY) 

C  INPUT 

C  VNIX  =  DERIVATIVES  OF  SHAPE  FUNCTION  W. R. T.  X,Y,Z 

C  OUTPUT 

C  VBE  =  MATRIX  B 

Q*********************************************************************** 


IMPLICIT  REAL*8  (A-H.O-Z) 
DIMENSION  V N I X ( 3 2 , 1) ,VBE(6,1) 
DO  10  1=1,6 


DO  10  J=1 . 96 
.10  VBE(I,J)=fl. 000 


3.5 


FORMATION  OF  MATRIX  B 


DO  20  1=1,32 
11=3*1-2 
12=11+1 
13=12+1 
KK=1 

Qi=y.Nix(i,^+] 


C2=VDIIXfI  (KK+l)) 
C3=VNIXf I  CKK+2)) 
VBE(lfIl)=Cl 


VBE(2  12  =C2 
VBE(3  I3}=C3 
VBEi  4  II  i=C2 
VBEi  4  12  =C1 
VBEi  5  I2)=C3 
VBE(5  I3)=C2 
VBE(6,Ilj=C3 
20  VBE(6  I3)=C1 


V  Ui-t  V  ,  4 

RETORT 
END 

SUBROUTINE  BTDB(VKE,VBF,VDE,IDLE, IMATD, NSYM) 

******  *************  *****&l+<i***i****l***1'*i**it1'4******it*i'****1'**1'*1'*it* 

C  TO  ADD  THE  PRODUCT  B(  i ). D. B  TO  THE  ELEMENT  MATRIX  K 
C  INPUT 

C  VBE  =  MATRIX  B 

C  VDE  =  MATRIX  D 

C  OUTPUT 

C  VKE  =  ELEMENT  MATRIX  K 

g*********************************************************************** 
DIMENS I0NRVKEf  1)(,  VBE(?Ma)d,  1)  ,VDE(  IMATD,  1)  ,T(576) 


10 

20 


DATA  ZERO/O.  0D07 
IJ=1 

IMAX=IDLE 
DO  40  J=1 . IDLE 
DO  20  I 1=1 , IMATD 
C=ZERQ 

DO  10  Jl=l, IMATD 


C=C+VDE(I1  J1)*VBE(J1,J) 
T(in=C 

F( NSYM. E Q v 0 )  IMAX=J 
1A> 


30 

40 


IF  ,  , 

DO  40  1=1 , IMAX 
C=ZERO 
DO  30  Jl=i, IMATD 
C=C+VBE(J1.I)*T(J1) 
VKEJIJj=VK£(lJ)+C 


RETURN 
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